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Abstract.  We  present  a  detailed  study  of  quantum  simulations  of  coupled  spin 
systems  in  surface-electrode  (SE)  ion-trap  arrays,  and  illustrate  our  findings 
with  a  proposed  implementation  of  the  hexagonal  Kitaev  model  (Kitaev  A  2006 
Ann.  Phys.  321  2).  The  effective  (pseudo)spin  interactions  making  up  such 
quantum  simulators  are  found  to  be  proportional  to  the  dipole-dipole  interaction 
between  the  trapped  ions,  and  are  mediated  by  motion  that  can  be  driven  by 
state-dependent  forces.  The  precise  forms  of  the  trapping  potentials  and  the 
interactions  are  derived  in  the  presence  of  an  SE  and  a  cover  electrode.  These 
results  are  the  starting  point  to  derive  an  optimized  SE  geometry  for  trapping 
ions  in  the  desired  honeycomb  lattice  of  Kitaev’s  model,  where  we  design  the 
dipole-dipole  interactions  in  a  way  that  allows  for  coupling  all  three  bond  types 
of  the  model  simultaneously,  without  the  need  for  time  discretization.  Finally,  we 
propose  a  simple  wire  structure  that  can  be  incorporated  into  a  microfabricated 
chip  to  generate  localized  state-dependent  forces  which  drive  the  couplings 
prescribed  by  this  particular  model;  such  a  wire  structure  should  be  adaptable 
to  many  other  situations. 


4  Author  to  whom  any  correspondence  should  be  addressed. 


New  Journal  of  Physics  13  (2011)  115011 
1 367-2630/1 1/11 501 1 


©  IOP  Publishing  Ltd  and  Deutsche  Physikalische  Gesellschaft 


2 


IOP  Institute  of  Physics  DEUTSCHE  PHYS1KALISCHE  GESELLSCHAET 


Contents 


1.  Introduction  2 

2.  Electrostatics  in  the  presence  of  conducting  planes  3 

2.1.  The  shielding  effect  of  the  cover  plane .  5 

2.2.  The  potential  due  to  the  surface  electrodes  (SEs) .  5 

2.3.  Dipole-dipole  interactions  between  trapped  ions .  7 

3.  Spin-spin  interactions  between  trapped  ions  8 

3.1.  Normal  modes  of  vibration .  8 

3.2.  State-dependent  forces .  9 

3.3.  Effective  spin-spin  interactions .  10 

4.  The  Kitaev  model  13 

4.1.  Model  and  implementation .  13 

4.2.  SE  trap  design .  16 

4.3.  Trap  depth  and  trap  loading .  17 

4.4.  Wires  for  magnetic  interaction .  18 

5.  Conclusions  21 

Acknowledgments  21 

Appendix.  Summary  of  the  used  coordinate  systems  21 

References  22 


1.  Introduction 

Ion-trap  systems  have  proven  to  perform  well  for  implementing  the  basic  elements  of  traditional 
quantum  computing,  where  evolution  is  described  in  terms  of  discrete  gate  operations, 
which  can  be  implemented  step  by  step  as  intermediate  states  are  irrelevant.  This  is  in 
contrast  to  quantum  simulations,  where  the  goal  is  to  simulate  the  continuous  evolution  of 
a  given  Hamiltonian.  While  the  initial  proposal  for  quantum  computing  with  trapped  ions 
relied  on  a  number  of  sequential  steps  to  mediate  effective  qubit  interactions  [1],  other 
approaches  [2-10]  achieve  interaction  between  the  internal  states  of  the  ions  via  constant 
Hamiltonians  and  therefore  allow  the  development  of  quantum  simulators  based  on  trapped 
ions  [7,  8,  11-15].  In  such  simulators,  interactions  between  trapped  ions  are  dominated  by  the 
Coulomb  potential.  For  this  interaction  to  affect  internal  states  (i.e.  the  qubits  or  pseudo-spins 
representing  the  effective  quantum  system  to  be  simulated),  state-dependent  forces  must  be 
applied  to  some  or  all  of  the  trapped  ions.  State-dependent  forces  can  be  achieved  through 
optical  ac-Stark  shifts  [5,  7,  15-19],  through  static  magnetic-field  gradients  in  combination 
with  homogeneous  radio-frequency  (rf)  fields  [4,  10,  14,  20]  or  with  rf  field  gradients  [21, 
22].  While  in  most  cases  the  Coulomb  interaction  is  considered  between  ions  in  a  self- 
assembled  single  chain  or  crystal,  coupling  of  independently  trapped  ions  has  recently  been 
demonstrated  [23,  24], 

For  quantum  simulations  with  ions  in  microtraps,  we  must  take  into  account  how  the 
presence  of  the  electrodes  modifies  the  Coulomb  interaction.  While  in  many  systems  this 
effect  is  negligible  (for  example,  in  the  surface  electrode  (SE)  setup  of  [23],  the  Coulomb 
coupling  was  found  to  be  enhanced  by  only  1.8%,  in  agreement  with  our  more  general  results 
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Figure  1.  Coulomb  interaction  between  two  charges  Q  and  Q'  (full  red  circles) 
in  the  presence  of  a  grounded  plane.  The  image  charges  (empty  red  circles) 
are  located  below  the  grounded  plane  and  carry  opposite  charge.  Interactions 
between  the  charges  (full  blue  arrow)  contribute  to  (1)  in  full,  whereas 
interactions  between  charges  and  image  charges  (dashed  blue  arrows)  contribute 
with  a  prefactor  of  1/2  [26]. 


in  section  2.3),  the  general  theory  developed  here  for  a  lattice  of  SE  microtraps  shows  that 
significant  modifications  to  free-space  couplings  are  possible.  Far  from  being  an  inconvenience, 
these  modified  interactions  can  be  used  to  design  quantum  simulations  with  specific  short-range 
effective  pseudo-spin  interactions,  which  we  illustrate  with  the  hexagonal  Kitaev  model  as  a 
concrete  example. 

The  remainder  of  this  paper  is  organized  as  follows.  In  section  2  we  present  a  Green’s 
function  approach  to  solving  electrostatic  problems  as  they  occur  for  SE  ion  traps  in  the  presence 
of  a  cover  electrode.  In  section  3,  we  derive  general  expressions  for  spin-spin  couplings  in  two- 
dimensional  (2D)  microtrap  arrays,  applicable,  for  example,  to  electric  coupling  to  light  fields 
or  magnetic  coupling  to  microwave  near- field  gradients.  In  section  4,  we  combine  all  these 
methods  to  show  how  the  hexagonal  Kitaev  model  [25]  can  be  implemented  with  an  array 
of  trapped  ions  on  an  optimized  SE  chip,  including  a  dedicated  wire  structure  that  could  be 
integrated  in  the  chip  to  simultaneously  mediate  the  couplings  along  three  distinct  bonds  by  the 
use  of  magnetic-field  gradients.  Finally,  the  appendix  gives  a  summary  of  the  used  coordinate 
systems. 


2.  Electrostatics  in  the  presence  of  conducting  planes 


The  electrostatic  interaction  between  charged  particles  close  to  conducting  surfaces  can  be 
strongly  modified  by  the  presence  of  the  conductors  [27].  In  the  idealized  geometry  of  a 
perfectly  conducting  grounded  electrode  plane  at  z  =  0,  the  total  electrostatic  energy  of  a  set 
of  charges  <2,  located  at  positions  r;  in  the  half-space  z-,  >  0  (see  figure  1)  is  [26] 


E 


c 

OO 


1 

4  7te0 


-E— + 

4-^  4zi 


Y  QiQjGoo(ri,rj ) 
‘<j 


(1) 


The  Coulomb  interaction  term  in  (1)  is  expressed  in  terms  of  the  Dirichlet  Green’s  function 
Goo,  which  can  be  found  from  the  free-space  Green’s  function  G(0)(r,  r ')  =  l/||r  —  r7||  by  the 
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Figure  2.  (a)  Sketch  of  an  SE  trap  with  a  grounded  cover  plane  positioned  at  a 
height  H  above  the  electrode  plane.  The  red  ring  electrodes  are  at  rf  potential, 
while  all  grey  areas  are  grounded.  For  static  interactions  or  interactions  varying 
slowly  compared  to  the  rf  period,  only  the  time-averaged  potential  contributes,  so 
for  our  purposes  the  situation  is  equivalent  to  two  completely  grounded  planes, 
(b)  The  interaction  energy  (3)  between  two  point  charges  at  same  height  h  over 
the  electrode  plane,  as  a  function  of  the  charge  separation  p,  in  the  presence  of 
a  cover  plane  at  height  H  =  100  h.  The  red  and  blue  parts  of  the  solid  curves 
are  computed  by  (4 a)  and  (4b),  respectively,  while  the  dashed  lines  illustrate  the 
approximate  behavior  given  by  (5). 


method  of  images  (see  figure  1), 


Goo(r ,  r') 


1 

yV  +  (z  ~  Z')2 


1 

vV  +  (z  +  z')2 


(2) 


where  p  =  sj(x  —  x’)1  +  (y  —  y')2  is  the  horizontal  distance  between  the  charges. 

In  the  following,  we  review  the  effects  of  a  grounded  cover  plane,  i.e.  a  second,  parallel 
conducting  plane  covering  the  electrode  plane  at  height  z  =  H  (see  figure  2(a)).  In  the  initial 
proposal  [28]  and  demonstration  [29]  of  SE  rf  traps,  the  conducting  surface  nearest  to  the 
trap  electrodes  was  theoretically  at  infinity  but  in  practice  a  part  of  the  surrounding  apparatus. 
It  has  been  suggested  that  adding  a  cover  plane  in  the  form  of  a  dc-biased  mesh  above  the 
electrodes  could  improve  trap  depth  [30].  In  addition  to  possible  benefits  of  providing  bias 
field  and  shielding,  the  cover  plane  could  have  more  practical  advantages,  namely  shielding  the 
trapping  region  from  fields  due  to  quasi-static  charges  on  insulators  in  the  vacuum  chamber  and 
establishing  a  more  well-defined  boundary  condition.  Further,  if  the  cover  plane  is  modified 
to  carry  rf  and  dc  electrodes  of  arbitrary  shape  in  the  same  way  as  the  electrode  plane,  the 
presented  formulae  can  be  used  directly  to  calculate  the  combined  electric  fields  generated  in 
this  ‘sandwich  trap’  geometry  (however,  if  optical  access  to  such  a  trap  geometry  is  achieved 
with  holes  and/or  fiber  optics  in  the  electrode  planes  [31],  the  present  full-plane  treatment  must 
be  adapted  [32]). 

Below,  we  first  modify  the  Green’s  function  (2)  to  include  the  cover  plane  and  illustrate 
that  a  cover  plane  at  height  H  leads  to  exponential  shielding  on  a  lateral  length  scale  of  H 
(section  2.1),  then  consider  its  effects  on  the  electric  potential  generated  by  SEs  (section  2.2) 
and  on  effective  dipole-dipole  interactions  between  vibrating  trapped  ions  (section  2.3). 
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2.1.  The  shielding  effect  of  the  cover  plane 


When  a  grounded  conducting  cover  plane  at  height  z  =  H  is  added  to  the  setup  of  figure  1,  the 
Coulomb  interaction  (1)  of  charges  located  between  these  two  planes  is  modified  to 


Ec 


1 

4jre0 


—  ^  Q*eH(Zi)  +  ^2  QiQjGH(ri,  r j) 
i  i  <j 


(3) 


Both  the  scaled  self-potential  eH(z)  and  the  Dirichlet  Green’s  function  GH  corresponding  to  the 
cover  plane  geometry  with  infinite  conducting  electrode  planes  at  z  —  0  and  z  =  H  can  be  found 
by  summing  over  an  infinite  sequence  of  mirror  planes;  and  in  the  absence  of  a  cover  plane 
( H  oo)  they  reduce  to  (1)  and  (2).  The  scaled  self-potential  is  eH(z)  =  —  [2y  +  \j/(z/ H)  + 
( 1  —  z/ H)]/ (AH)  —  E  +  0(z2/H3)  in  terms  of  Euler’s  constant  y  =0.577  216. . .  and  the 
digamma  function  iff  (a)  =  T'(a)/  T(a).  The  Dirichlet  Green’s  function  is 


OO 

GH(r,r')  =  Y,  G00(r+2/x//z,r/) 

fl=  —  00 


(4  a) 


m 


where  Kq  is  the  modified  Bessel  function  of  the  second  kind.  The  second  form  (Ah)  is  obtained 
by  solving  the  Laplace  equation  in  cylindrical  coordinates  [27].  Both  forms  converge  for  all 
parameters  (p,  z,  z'),  but  while  (Aa)  converges  faster  when  ||r  —r' ||  <  H,  (Ab)  is  more  suitable 
if  ||  r  —  r'  ||  >  H ,  in  particular  for  p  H,  as  discussed  below. 

The  Coulomb  interaction  energy  GH(p,  z,  z')QQ' /(Attcq)  between  two  charged  particles 
in  an  SE  trap  depends  on  the  horizontal  separation  p,  as  illustrated  in  figure  2(b).  To  illustrate 
this  interaction  energy  we  take  the  particles  to  be  at  the  same  height  h  above  the  electrode  plane, 
and  the  cover  plane  height  H  to  be  much  larger  than  h.  When  p  <<P  H .  we  expect  the  cover 
plane  to  be  irrelevant,  so  that  the  interaction  is  described  by  a  single  image  charge:  it  falls  off 
as  p~l  when  p  <  h  (where  the  electrode  plane  is  irrelevant)  and  as  p~3  thereafter,  as  described 
by  (2).  When  p  >  H  the  cover  plane  becomes  important  and  the  asymptotically  dominant  form 
is  the  first  term  of  the  resummation  (Ab),  so  that  the  presence  of  the  cover  plane  leads  to  an 
exponential  shielding  at  the  length  scale  of  the  cover  plane  height,  as  illustrated  in  figure  2. 
Summarizing, 


GH(p,  h,  h) 


1  fp, 

2  h2/p\ 


e-np/H 


for  p  <^h, 
for  h  <£  p  «  H, 

for  pyp  H. 


(5) 


2.2.  The  potential  due  to  the  surface  electrodes  (SEs) 

The  contribution  to  the  total  potential  from  the  structured  electrodes  in  the  z  =  0  plane  can  be 
computed  as  an  integral  over  the  electrode  plane: 

&(r)=  f  G{^\r ,  r')<J>(r')djc'dy/,  (6) 

Jz'= o 
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where  we  have  introduced  a  ‘surface  Green’s  function’  G^\r,  r')  =  (r,  /•/)|_,_0-  In  the 

absence  of  a  cover  plane,  the  surface  Green’s  function  was  found  to  be  [32] 

Gg>(r,r')  =  GgVZ)=  Z  (7) 

2tt  ( p2  +  z 2) 

with  the  geometric  interpretation  that  the  potential  at  r  due  to  an  electrode  at  potential  4>0 
is  <f>0/2 7T  times  the  solid  angle  spanned  by  the  electrode  as  seen  from  r  [33].  Alternatively, 
the  electric  field  at  r  is  proportional  to  the  magnetic  field  that  would  be  observed  if  a  current 
were  flowing  along  the  edge  of  the  electrode  [33,  34].  For  electrode  configurations  that  are 
translationally  invariant  in  the  x -direction,  the  system  can  be  described  by  conformally  mapping 
the  upper-half  yz-plane  (z  >  0)  to  a  disc  [35].  Analogous  to  (4 a)  and  (4b),  we  have  for  the 
general  case  including  a  cover  plane, 

OO 

Gf(p,z)=  J2  G™(p,z  +  2liH)  (8 a) 

fl=  —  00 


.  OO 

=  °^z)-—2  £  o'+ixo'+2)(^)G,Q, 

.7=1,3,... 


(8c) 


where  £  is  the  Riemann  zeta  function,  Pj  are  Legendre  polynomials  and  s=  ||r  —  r'||  = 
y/p2  +z2.  Forms  (8a)  and  (8b)  converge  for  all  (p,  z),  but  while  (8a)  converges  faster  when 
s  <  H ,  (8b)  is  more  suitable  if  s  >  H .  Form  (8c)  is  restricted  to  s  <  2 H  and  is  most  useful  for 
s  <$C  H.  Similar  to  (5)  we  find  the  approximate  behavior 


G(h\P,z) 


8  ttH2 

1 


sin  (f ) 


,-np/H 


for  p  <£  z, 
for  p  »  H, 


(9) 


where  \lr'(a)  =  T"(a)/T(a) —  ^2(a)  is  the  first  derivative  of  the  digamma  function.  We 
conclude  that  the  influence  of  any  SE  is  exponentially  damped  at  distances  larger  than  H,  which 
is  advantageous  for  the  experimental  construction  of  quasi-infinite  surface  microtrap  lattices  in 
that  it  reduces  the  influence  of  the  inevitable  electrode  boundary:  at  any  point  further  than  H 
away  from  the  edge  of  the  electrode  and  cover  plane,  the  trap  will  look  as  if  the  electrode  were 
infinitely  large. 

Since  the  surface  Green’s  function  only  depends  on  the  x-  and  y -coordinates  through 
r  —  r' ,  (6)  is  a  folding  integral  (convolution)  [32]  and  can  be  rewritten  as  a  product  of  the 
Fourier- transformed  quantities,  (kx ,  ky,z)  =  G(^](kx,  ky,  z)&(kx,  ky,  0),  with 

1  f00 

®(kx,  ky,  z)  =  —  ®(x,y,  z)e-*k'x+k>y) dx  dy  (10) 

j _00 


and  a  similar  expression  for  the  Fourier-transformed  Green’s  function.  The  latter  is  cylindrically 
symmetric  (k  =  k2  +  k2), 


Gf(k,z) 


sinh(k//  —  kz) 
sinh(k//) 


for  H  — >  oo, 


(11) 
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and  allows  a  rather  intuitive  interpretation.  All  solutions  of  the  Laplace  equation  with 
the  horizontal  wavevector  {kx,ky}  are  of  the  form  c,ikxX+k'y>  (a+c+kr'  +  <x_erkz)\  the  Green’s 
function  (11)  gives  the  unique  solution  which  has  unit  amplitude  on  the  electrode  plane 
[G(jj\k,  0)  =  1]  and  zero  amplitude  on  the  cover  plane  [  (k,  H)  =  ()].  Therefore,  (11)  gives 

the  unique  extension  of  a  unit- amplitude  potential  plane  wave  from  the  z  =  0  plane  into  the 
z  >  0  half-space  which  satisfies  the  boundary  condition  of  vanishing  amplitude  on  the  cover 
plane.  The  fact  that  the  momentum-space  representation  of  the  surface  Green’s  function  (11) 
can  be  written  without  infinite  sums  greatly  simplifies  the  description  of  infinite  lattices  of  SE 
microtraps  [36]. 


2.3.  Dipole-dipole  interactions  between  trapped  ions 


Trapped-ion  quantum  simulators  couple  internal  degrees  of  freedom  of  the  ions  (typically 
hyperfine  states  or  metastable  D-states)  through  a  state-dependent  coupling  to  shared  vibrational 
degrees  of  freedom  [1,  6,  8,  14,  20,  21]  (see  section  3).  A  crucial  ingredient  of  these  couplings 
is  the  precise  nature  of  the  Coulomb  interactions  between  the  ions.  Here  we  address  the  details 
of  this  latter  point,  since  it  will  determine  how  to  construct  a  quantum  simulator  of  a  desired 
system,  as  exemplified  in  section  4. 

We  consider  the  regime  of  ‘stiff’  ion  trapping  [6],  where  the  Coulomb  interaction  is 
relatively  small  compared  to  the  trapping  potential,  and  we  can  interpret  the  normal-mode 
dynamics  of  the  ion  crystal  as  that  of  a  set  of  local  harmonic  oscillators  that  are  weakly  coupled. 
The  ion  trapping  potential  defines  a  set  of  local  eigenmodes  for  the  /  th  ion  corresponding 
to  vibration  in  three  orthogonal  directions  m  ‘  (with  \\mf  ||  =  1  for  p.  =  1,  2,  3)  around  an 
equilibrium  position  R0i.  In  what  follows,  we  use  these  directions  to  parameterize  the  position 
of  the  f  th  ion  as 

3 

A  =  R()i  +  ^  rim'..  (12) 

ii=i 


The  total  Coulomb  energy  of  a  set  of  N  charges  is  given  in  (3),  and  the  leading-order  term 
that  couples  the  motion  of  the  ions  is 


j-i  coupling 


l 

4ttc0 


N  3 

J2  J2  QiQjrfr'jm?  -  Vy jGniRoi,  R0j)  -mvj. 

i<j  ii. v=l 


(13) 


Since  we  are  mainly  interested  in  near(est)-neighbor  interactions,  we  evaluate  this  expression  in 
terms  of  the  infinite  sum  over  image  charge  pairs  (4a),  rather  than  the  resummed  form  (4b): 

OO 

m  ■  VV'G//(r,  r')  •  m  =  ™  ■  VV'G00(r  +  2pHz.  r')  •  in’,  (14) 

/x=— OO 


where  the  explicit  dipole-dipole  coupling  is  given  by  the  expression  without  a  cover  plane, 


m  ■  VV'Goolr,  r’)  ■  mf 


m  ■  m!  —  3(m  ■  n)(mr  ■  n)  m  ■  m’  —  3(m  ■  n)(m’  •  n) 


-  r'l|3 


-  r’  1 1 3 


(15) 


in  terms  of  it  =  (r  —  r')/||r  —  r'||,  n  =  (r  —  r')/||r  —  r'||,  and  the  mirrored  quantities  r '  = 
r’  —  2 (rr  ■  z)z  and  in'  =  in'  —  2(m'  ■  z)Z-  The  first  term  of  (15)  is  the  well-known  dipole-dipole 
interaction,  while  the  second  term  is  the  correction  due  to  image  charges  in  the  electrode. 
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In  order  to  illustrate  the  behavior  of  the  dipolar  interaction  (15)  in  the  close  proximity  of 
a  conducting  electrode  plane,  we  again  consider  two  ions  located  at  equal  height  h  above  the 
electrode  plane,  spaced  by  a  distance  p  along  the  x-axis  and  in  the  absence  of  a  cover  plane. 
If  we  assume  that  both  ions  vibrate  along  axes  m  =  mr  that  are  parallel  to  the  laboratory-frame 
coordinate  axes,  then  we  find  that  the  presence  of  the  electrode  plane  can  either  increase  or 
decrease  the  dipolar  coupling  strength: 


m  =  m'  =  x  : 


m  =  m’  =  y  : 


m  —  m'  —  z  : 


x-WG0o(hz,hz  +  px)-x  = - - 


p3(2h2  —  p 2)  1  p»2/?  24  h2 


1  + 


p3  (p2  +4h2)5/2  J  p 


5  ’ 


1 


y  ■  VV'GqoC/iz,  hz  +  px)  ■  y  =  +— 


P 


p3  L  (p2  +  4fi2)3/2_ 


p»2/i  6  h2 


+- 


P 


5  ’ 


z  ■  V  V'Goo(fi£,  hz  +  px)  ■  z  =  +^r 

p3 


1- 


p3(Sh2  —  p2)1 
(p2  +  4/i2)5/2 


p»2/i  2 


(16) 


Thus  we  see  that  by  choosing  the  directions  of  vibration  m  in  particular  ways,  we  can  use  the 
presence  of  the  electrode  plane  to  make  the  dipolar  interactions  fall  off  with  the  fifth  power  of 
distance  instead  of  with  the  third  power  (which  is  the  case  in  the  absence  of  any  conducting 
planes),  as  long  as  the  ion  oscillation  frequency  is  low  enough  to  avoid  the  effects  of  retardation 
and  dissipation.  The  relevant  length  scale  that  determines  whether  or  not  the  electrode  plane 
has  a  strong  influence  on  the  dipole-dipole  coupling  is  p  ~  2h,  similar  to  figure  2;  for  even 
farther  separations  (p  >  H),  we  find  exponentially  damped  dipole-dipole  couplings  due  to 
the  shielding  effect  of  the  cover  plane  (see  section  2.1).  These  rapid  dampings  can  be  used 
to  construct  lattice  simulation  models  with  nearly  local  interactions,  which  is  a  desirable  feature 
since  many  spin  models  from  condensed-matter  physics  are  formulated  in  terms  of  such  local 
(e.g.  nearest-neighbor)  couplings. 


3.  Spin-spin  interactions  between  trapped  ions 

This  section  derives  how  state-dependent  forces  can  induce  pseudo-spin  interactions  between 
neighboring  ions  through  the  Coulomb  potential.  While  this  effect  is  well  known  in  principle  [1], 
we  show  how  these  effective  interactions  are  constructed  in  a  lattice  of  ions  without  the 
need  for  time-slicing  (‘Trotterization’  [37]).  Further  we  show  that,  to  lowest  order,  the 
effective  interaction  strengths  are  proportional  to  the  real-space  Coulomb  coupling  strengths, 
an  observation  that  greatly  simplifies  the  design  of  lattice-based  quantum  simulators  (see 
section  4). 

3.1.  Normal  modes  of  vibration 

For  small  oscillation  amplitudes  r!‘  the  coupled  harmonic  motion  of  N  ions  in  a  lattice  can  be 
described  by  considering  the  local  trapping  potential  curvatures  (the  second  derivatives  of  the 
ion  trapping  pseudo-potential  with  respect  to  position)  around  the  equilibrium  positions  Roi  and 
including  the  Coulomb  couplings  between  ions  in  separate  wells  to  second  order  [38].  If  we 
assume  that  all  excursions  r,  —  R0i  are  already  written  in  terms  of  the  ‘bare’  eigenmodes  of 
the  isolated  local  trapping  potentials  (with  ‘bare’  frequencies  wi/x),  as  in  (12),  then  the  potential 
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energy  of  the  ions  of  mass  M  is 


V  = 


N 


N 


£E“?X>2+£  £ 


/XV  /X  V 

y.  r  r  ■ 
•  n  11 


i = 1  i±=  1 


*',7=1  M,v=l 


(17) 


where  y£v  =  0  and  ■  V ,  V  jGH(R0i,  Roj)  ■  my  see  (13).  This  quadratic  potential 

energy  can  be  diagonalized  using  coefficients  Oitim  such  that  the  real-space  displacements  can 


be  written  as 

3  N 

3  N 

N  3 

t"/  —  ^  '  f//  lull  q rn , 

with  ^  ^  OifimO jvm  —  SijSf^y  and 

^  ^  ^  ^  0 i  fim  0 i  fim'  —  &mm'  i 

(18) 

m=l 

- — .  a  at  i 

m  —  1 

i=l  f±=  1 

and  V  =  Y^m=\  ’n  terms  of  the  lattice  normal-mode  amplitudes  qm  and  their 

frequencies  com.  We  quantize  these  normal  modes  through  qm  i->-  qm  =  qom  (am  +  am)  with  q0m  = 


y/fi/ilMcOm)  and  the  usual  commutation  relations  [am ,  am,]  =  8mm>.  We  will  work  in  the  ‘stiff’ 
lattice  limit,  where  we  assume  that  the  ‘bare’  trap  frequencies  Vi  are  equal  for  all 

ions5  and  the  Coulomb  couplings  between  ions  do  not  significantly  mix  local  modes  with 
different  values  of  /z.  This  is  the  case  when  (i)  the  bare  trap  frequencies  are  sufficiently  far 
apart:  \y^v\  <?C  [dr,  —  drv\  V  i,  j,  /x^y,  and  (ii)  the  vibrational  bands  are  sufficiently  narrow: 
\y^\  «:  miny  or  —  oo2v\  V  i,  j,  //.  In  this  limit,  the  normal  modes  of  the  lattice  separate  into 
three  disjoint  sets  indexed  by  /z  e  {1,  2,  3},  each  containing  N  normal  modes  with  frequencies 
close  to  the  corresponding  co^. 


3.2.  State-dependent  forces 

In  addition  to  the  Coulomb  couplings,  which  are  always  on  and  define  the  coupled  vibrational 
eigenmodes  of  the  trapped  ions,  we  can  experimentally  introduce  fields  that  couple  to  internal 
states  of  the  ions.  Examples  of  such  interactions  are  electric  or  magnetic  dipole  couplings, 
Raman  couplings  or  electric  quadrupole  couplings  to  laser  or  microwave  fields.  In  the  following 
general  treatment,  we  assume  that  there  is  a  coupling  between  (a)  classical  external  field(s) 
effectively  oscillating  with  angular  frequency  coi,  and  two  internal  states  of  each  ion,  forming 
an  effective  two-level  (spin- 1/2  or  pseudo-spin- 1/2)  system.  Irrespective  of  the  type  of  induced 
coupling,  the  coupling  operator  of  the  ith  ion  in  its  (pseudo)spin-l/2  subspace  can  be  expressed 
as  a  linear  combination  of  the  identity  operator  a0(,)  and  the  Pauli  matrices  by(},  l  e  {X,  Y,  Z}, 
expressed  in  a  quantization  coordinate  frame  whose  axes  are  given  by  the  orthonormal  vectors 
X,  Y,  Z  (see  the  appendix).  The  coupling  Hamiltonian  can  thus  be  very  generally  expressed  as 

N 

T-L\  ^  ^  ^  [cf  cos(<z>p  +</>c°)  +  (ji  -  R0i )  •  4°  cos(tuif  +  0s(())]ct£(O,  (19) 

i= 1  Xe{0,  X.Y.Zj 

where  we  have  performed  a  first-order  expansion  in  the  ion  positions  assuming  small  oscillation 
amplitudes.  Any  type  of  spin- 1/2  coupling  that  is  used  with  trapped  ions  (including  effective 
couplings  to  pseudo-spin  degrees  of  freedom)  can  be  brought  into  this  form,  where  terms 

5  This  equality  can  be  relaxed  to  the  condition  that  for  the  modes  that  are  used  for  inducing  spin-spin  couplings, 
the  spread  of  the  bare  frequencies  is  much  smaller  than  the  dominant  Coulomb  couplings. 
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with  non-vanishing  prefactors  cp  and  sy  are  referred  to  as  ‘carrier’  and  ‘sideband’  terms, 
respectively.  The  phases  can  absorb  differences  in  the  details  of  driving  fields:  while  stationary 
fields  in  general  have  ffp—ffp,  travelling  waves  (e.g.  light  fields)  are  characterized  by 

4>f  =  08°  ±  I' 

As  an  example,  the  coupling  of  physical  spins  to  a  magnetic  field  is  found  by  expanding 
their  magnetic  dipole  operators  as  ft{,)  =  g(l)  pR(dP  X  +  dp}Y  +  opZ)  in  terms  of  the  Bohr 
magneton  /iB  and  the  g-factors  g<0;  for  small  ion  excursions  the  coupling  Hamiltonian  to  the 
magnetic  field  EL\  =  —  JPj  fi(,)  ■  B(r  ,  )  cos  (cop  +  0)  can  thus  be  expressed  in  the  form  of  (19) 


with 

c(0  _ 

f°- 

for  l  =  0, 

ce  ~ 

-g^pBX  ■  B(R0i), 

for  l  =  X, 

and  similarly  for  £  =  Y,  Z, 

A0  _ 

o, 

for  l  =  0, 

si  — 

-g«ffBV[X-B(R0i)l 

for  £  =  X, 

and  similarly  for  £  =  Y,  Z 

and  cj)p  =  4>P  =  cf)  V  i .  We  stress,  however,  that  the  above  form  of  the  magnetic  dipole  operator 
does  not  apply  to  pseudo- spins  for  their  effective  interactions  with  external  fields;  see  section  4.4 
for  an  example  involving  pseudo-spins. 


3.3.  Effective  spin-spin  interactions 


Inserting  the  lattice  normal-mode  expansion  (18)  and  (12)  into  (19),  we  can  write  the  interaction 
Hamiltonian  as 

Nr  3  N 


*■  =  £  E 

i= 1  ee{0 ,X,Y.Zj  u 


cf  COS  (cop  +  cj>P)  +  YP  2  tiQ.imt(dm  +  a]n)  cos  (cop  +  0s(l)) 


m= 1 


(0 


(21) 


where  we  have  dropped  the  approximation  symbol  and  introduced  £2 


iml 


sy .  It  is  common  to  transform  into  the  interaction  picture  to  assess  the  dynamics  induced 
by  such  an  interaction  Hamiltonian.  In  this  transformation,  the  field-free  Hamiltonian 
'Ho  =  Yl3m=i^ia>m(alldm  +  p)  +  h(o^iJ2f=1&P  leads  to  a  time  dependence  of  the  operators 
in  (21): 


—  W3  O  nff 

2ft  3—!iL= 1  '-'ifj-m"1, 


dm  i— >■  dmc  djn  dllc]0>mt, 

^  °x  cos(mn0  -  dp  sin(cun0, 


o. 


(0 


!->•  O. 


(0 


er 


(0 


!->•  O 


-  (0 


a 


(0 


i->-  a 


-(«) 


sin(cu|4,f)  +  dy }  cos (co^t). 


(22) 


The  terms  involving  d{p  and  dp  can  lead  either  to  spin  flips  without  affecting  the  motion 
(‘carrier’  transitions,  mediated  by  c(pd(p  and  cpdp  and  resonant  at  the  frequency  difference 
a>ji  between  the  pseudo-spin  states)  or  to  interactions  that  couple  spins  and  motion  (‘sideband’ 
transitions  or  Mpl mcr-Sprcnscn  interactions  [3],  mediated  by  sp dp  and  s(pdP  and  resonant 
around  &qi±<uM);  the  latter  will  dominate  if  they  are  not  driven  too  strongly  and  \o>\  — 
co^  ±  d>P  <$C  \cjO\  —  co^\  for  one  of  the  signs  in  ±.  Here,  we  concentrate  instead  on  a  drive 
with  frequency  \a>i  —  o)3[  <$C  fo>3  <$C  <«t;,  close  to  one  of  the  three  bare  eigenfrequencies  of  the 
uncoupled  ion  sites  (we  have  chosen  p  =  3  without  restricting  generality).  In  this  case,  we  can 
neglect  the  terms  in  sp  and  .sy’  as  they  are  far  off-resonant,  and  all  c(p  by  careful  design  of  the 
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experiment  (see  section  4.4).  The  interaction  Hamiltonian  in  the  interaction  picture  thus  reduces 
to  a  coherent  drive 


U |nt 


N  3  N 

fi  J2  E(«me-i(J"'-^  +al^Smt-^ 

i=  1  m  —  1 


)(Qim0a^  +  Qimz^‘)) 


OK 


(23) 


after  a  second  rotating- wave  approximation,  with  the  detunings  Sm  =  wm  —  oj\.  Equation  (23) 
can  be  exactly  integrated  via  a  Magnus  expansion  [39, 40]  to  yield  the  unitary  evolution  operator 


U{m(t)  =  exp 


r  N  3N  i 

EE1'6 

*"  1  =  1  m= 1 


t 


-e  '4,!-'>all(Qim0a^)  +  QimZa^)  -  h.c. 


exp 


r  N 

i  J2  cos  0s; 

L  i,j= 1 


3N  -on 

v  Y^/n  ~  ( 0  I  o  iWun  ~  O')  ,  o  ~ O' ) \  ^  sin  (<5m  / ) 

X  /  T^^'’mQCr0  '  )\»£jm0&0  '  ^^jinZ^z  /  ^9  > 

m=l  m 

(24) 


with  =  0(0  _  00') _  The  first  exponent  describes  a  set  of  time-dependent  coherent 
displacements  to  all  normal  modes  that  can  entangle  the  motion  with  the  internal  pseudo¬ 
spin  states.  The  second  exponent  constitutes  a  phase  that  depends  on  pairs  of  spins  and  can  be 
interpreted  as  a  spin-spin  interaction.  For  a  faithful  simulation  of  interacting  spins  it  is  desirable 
that  (A)  the  first  term  should  be  as  close  as  possible  to  the  identity  operator,  in  order  to  avoid 
populating  vibrational  excitations,  and  (B)  the  second  term  should  provide  sizable  phases  for 
desired  inter-ion  couplings  with  i  ^  j,  as  these  represent  the  spin-spin  interactions.  It  can  be 
shown  from  the  expression  above  or  by  the  use  of  a  canonical  transformation  [6]  that  (A)  can 
be  approximately  met  as  long  as  <$C  |<5m|  for  all  (i,  m,  £). 

The  above  restrictions  do  not  limit  the  time  scale  for  simulations,  as  long  as  one  assumes 
that  sufficiently  strong  couplings  can  be  induced  by  lasers  or  microwave  field  gradients. 
However,  the  energy  scale  of  nearest-neighbor  Coulomb  interactions  also  plays  an  important 
role  in  determining  simulation  time  scales,  but  this  dependence  is  hidden  in  the  normal-mode 
coefficients  Oilim  of  (18).  To  illustrate  this  point,  we  assume  that  £2/m0  =  0  for  all  normal  modes 
m  and  sites  i  (see  (20)  for  an  example),  but  what  we  show  below  also  holds  for  more  general 
cases.  Assuming  negligible  displacements  (see  point  (A)  above)  the  unitary  evolution  operator 
thus  simplifies  to 


U\{t)  =  exp 


N 


3  N 


En' 


imZ^jmZ  ^ 


8mt 


sin  (8mt) 


ij= 1 


m= 1 


8i 


for  t  »  sup  \8m]\.  (25) 

m 

In  the  ‘stiff’  lattice  limit  (see  section  3.1),  we  choose  the  drive  frequency  <x>i  close  to  one  of  the 
bare  frequencies,  say  op ,  such  that  the  detunings  Sm  will  be  much  smaller  for  normal  modes 
in  this  set  than  for  the  other  normal  modes;  consequently,  the  sum  over  modes  m  in  (25)  can 
be  restricted  to  an  ‘active’  set  of  N  modes  clustered  around  0)3.  If  we  further  choose  the  drive 


:  exp 


N 

kE 

U=1 


a j  ]  cos  (p'J 


3  N 

E 

m= 1 


jmZ 
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frequency  such  that  <$3  =  (03  —  a>i  is  much  larger  than  the  spread  of  the  normal  mode  frequencies 
in  the  active  group,  the  series  expansion 


1  =  1  3 


+  0[«-a>32)2] 


63  20)363 


together  with  the  relations  in  the  ‘active’  group  of  normal  modes 


N 


^  '  Oj^mO j^m  —  Sjj , 


m= 1 


N 

E 

m= 1 


Oi3mO  j3m(^m  ^3  )  —  Y( 


33 

ij 


(26) 


(27a) 


(21b) 


simplifies  the  unitary  evolution  (25)  to 
u/oV 


U\(t)  exp 


4  h% 


E<m.!'sz)2 


i  =  l 


x  exp 


i<7oV 


N 

8^2o)36? 

j  1,7=1 


X 


33  cos  4>lJ  (m] 


s{z)(m 


_0')\  -  (0  -  O') 
z  ;<rz  ^2 


(28) 


where  we  have  further  approximated  qom  ~  g03  =  s/ti/ (2M 0)^).  The  first  term  in  (28)  is  a  global 
phase;  it  is  the  second  term  that  mediates  an  effective  spin-spin  coupling  on  the  lattice  of  ions. 
It  can  be  interpreted  as  the  evolution  under  an  effective  spin-spin  coupling  Hamiltonian  due  to 
the  driving  of  local  mode  /x  =  3, 


U 


eff 
33  ' 


N 

Ex 

t,j= 1 


33  -  (0-0') 


o7  o 


(29) 


where  the  effective  spin-spin  coupling  coefficients  are 


qliYi* cos  €J  (mi  ■  sz)(ffj]  •  s'z’) 
Sfid>j,S^ 


(«h 


,0'h 


j  33  _ 

Jij  ~ 


(30) 


We  conclude  that  to  lowest  order  the  strength  of  the  spin-spin  coupling  between  two  ions  is 
determined  by  the  geometric  overlaps  ( m 3  ■  s(%)  and  (/n3  •  s 1 )  of  the  ‘active’  local  modes  of 
vibration  with  the  direction  of  the  state-dependent  force,  as  well  as  by  the  real-space  Coulomb 
coupling  strength  y;33  between  the  ions  moving  along  these  local  modes  (see  (13)  and  (17)). 
The  relative  phases  4>'J  of  the  driving  forces  can  be  used  to  modulate  the  coupling  strengths. 
These  observations  are  used  in  section  4  to  construct  a  quantum  simulator  on  a  lattice  of  trapped 
ions.  Equations  (29)  and  (30)  faithfully  describe  the  evolution  of  the  system  under  the  following 
conditions: 


•  The  vibrational  band  structure  of  the  trapped  ions  must  consist  of  clearly  distinct  bands 
which  can  be  addressed  individually  (see  section  3.1  for  the  conditions  for  ‘stiff’  trapping). 

•  The  state-dependent  force  must  be  driven  at  a  small  detuning  8  from  one  of  these  bands:  8 
must  be  large  enough  such  that  (26)  is  valid  for  this  band,  but  small  enough  such  that  the 
contributions  to  (30)  from  other  bands  are  negligible. 
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•  The  amplitude  of  the  state-dependent  force  must  be  small  enough  such  that  it  does  not 
significantly  excite  ion  vibrations.  The  condition  |  £2unt  |  <$C  \Sm  |  given  above  implies  that  the 
state-dependent  forces  must  be  weaker  than  the  force  scale  F  ~  fi\S\/q0  for  the  addressed 
band. 

The  above  arguments  can  be  made  in  a  very  similar  fashion  for  a^]a(xj)  and  dy'dy' 
interactions  by  considering  the  interaction  Hamiltonian  (21)  with  u>\  ~  ±  cu/x  in  the 

appropriate  basis  |±)  =  (|f)  ± elx ||))/\/2,  where  the  considered  spin-spin  interaction  is 
diagonal.  The  only  slight  complication  can  arise  from  carrier  terms  proportional  to  c'^  and  cy’ 
that  are  detuned  by  roughly  the  motional  eigenfrequencies.  For  detunings  from  the  sidebands  on 
the  order  of  the  dipole-dipole  interactions  and  correspondingly  small  drive  strengths,  however, 
these  carrier  terms  can  be  safely  neglected. 

To  summarize  this  section,  we  have  considered  general  effective  spin-spin  interactions  in 
the  limit  of  ‘stiff’  ion  trapping.  We  have  shown  that  even  in  a  lattice,  the  spin-spin  coupling 
strength  of  any  two  ions  depends  on  the  dipolar  Coulomb  coupling  between  these  two  ions.  To 
avoid  appreciable  entanglement  between  (pseudo)spins  and  ion  motion,  the  detunings  of  driving 
fields  need  to  be  larger  than  the  couplings  they  induce.  This  latter  finding  agrees  with  other  work 
on  simulation  with  trapped  ions  [8].  At  the  same  time,  the  detunings  cannot  be  much  larger  than 
the  couplings  between  nearest  neighbors,  which  determine  the  finer  structure  of  the  normal¬ 
mode  spectrum  around  the  frequencies  of  the  uncoupled  (‘bare’)  motion  of  an  ion  tightly  bound 
in  one  of  the  trapping  wells  (see  section  4  for  a  concrete  example).  These  requirements  impose 
stringent  bounds  on  the  time  scales  necessary  to  perform  simulations.  For  example,  in  [23] 
two  ions  at  a  distance  of  40  jim  exhibited  an  exchange  splitting  of  approximately  3  kHz,  barely 
sufficient  to  demonstrate  a  few  energy  exchanges  before  ion  heating  profoundly  altered  the 
motion.  Simulations  that  need  to  progress  adiabatically  with  respect  to  this  exchange  period  will 
therefore  be  experimentally  challenging  and  may  require  reducing  anomalous  heating  below 
what  was  measured  in  [23]. 

4.  The  Kitaev  model 

As  an  example  of  how  to  use  the  results  of  sections  2  and  3  in  the  design  of  a  quantum  simulator, 
we  construct  an  implementation  of  the  hexagonal  Kitaev  model  [25]  with  microtrapped  ions.  In 
its  ideal  form,  this  exactly  solvable  2D  spin  model  has  a  topologically  ordered  ground  state  with 
anyonic  excitations,  which  makes  it  extraordinarily  interesting  for  study  in  a  quantum  simulator 
with  individual  access  to  the  constituent  degrees  of  freedom. 

4.1.  Model  and  implementation 

The  hexagonal  Kitaev  model  [25]  has  the  Hamiltonian 

3^ Kitaev  — 1  Jx  ^  °X°X  >  *  J¥  ^  ^Y^y'  ~  JZ  ^  (31) 

X— links  Y— links  Z— links 

defined  on  a  honeycomb  lattice  of  spin- 1/2  particles,  where  the  laboratory-frame  bond  vectors 
refer  to  figure  3: 

Ax  =  d{0,  1,0],  Av  =  d{V3,  —1, 0}/2,  Az  =  d{- V3, -1, 0}/2.  (32) 
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Figure  3.  Dipole-dipole  interactions  (15)  with  the  central  (green)  site  due  to  the 
vibrational  directions  of  (33)  for  the  case  of  /x  =  X,  expressed  as  percentages 
of  the  dominant  coupling  (equal  to  1.99  x  Q2/(4ne0d3),  see  (34);  values  below 
1%  of  this  are  not  shown).  Figures  4  and  5  show  the  vibrational  band  structure 
induced  by  these  couplings.  The  couplings  for  Y  (Z)  are  found  by  rotating  this 
figure  by  120  (240  )  clockwise.  Red  and  blue  wires  are  described  in  section  4.4. 


In  this  way  the  Hamiltonian  (31)  associates  each  real-space  bond  direction  (32)  with  a  spin 
quantization  direction;  however,  it  is  important  to  keep  in  mind  that  the  bond  directions  and 
the  associated  spin  quantization  directions  are  not  a  priori  related  (see  appendix).  The  ions  are 
located  on  two  sublattices  L0  and  L.,  as  shown  in  figure  3.  Neighboring  ions  are  a  distance  d 
apart. 

The  form  of  (31)  is  exactly  that  of  (29)  summed  over  three  concurrent  driving  force  fields. 
As  these  driving  fields  will  be  at  very  different  frequencies,  they  can  be  applied  simultaneously 
in  order  to  drive  the  full  Hamiltonian  (31).  What  is  therefore  needed  in  order  to  implement  the 
Kitaev  model  is  a  set  of  ‘bare’  vibrational  directions  of  the  ions  such  that  the  couplings  y^ v, 
and  therefore  the  effective  spin-spin  couplings  (30),  match  the  particular  geometry  of  the  three 
terms  in  (31). 

We  choose  the  ion  trapping  height  to  be  half  of  the  inter-ion  distance,  h  =  d/2,  and  the 
orthonormal  principal  axes  of  vibration  for  ions  on  the  two  sublattices  as 

mx  =  {0,  2,  V2}/V6,  ml  =  {0,  -2,  V2}/V6, 

mY0  =  {v/3,  -1,  V2}/V6,  mY,  =  {-V3,  1,  V2}/V6,  (33) 

ml  =  {-V 3,-1,  V2J/V6,  ml  =  {V 3,  1,  V2J/V6. 

This  particular  choice  of  axes  of  vibration  has  the  property  that  dipole-dipole  couplings  of 
the  sort  of  y (i.e.  coupling  the  mf  vibration  of  the  ion  at  R0i  with  the  m1'  vibration  of  the 
ion  at  R()j)  are  strongly  dominated  by  the  nearest-neighbor  couplings  required  by  the  Kitaev 
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Figure  4.  Vibrational  band  structure  due  to  the  dipole-dipole  interactions 
of  vibrations  along  one  of  the  sets  of  axes  in  (33)  (/x  =  X,  see  figure  3), 
corresponding  to  the  density  of  states  shown  in  figure  5.  In  the  lower  band  (left) 
neighboring  ions  move  out  of  phase;  in  the  upper  band  (right)  they  move  in 
phase.  The  first  Brillouin  zone  is  drawn  in  black.  Frequencies  (colors)  are  given 
in  units  of  co0  (see  figure  5). 


model  (31),  shown  in  figure  3  for  /x  =  X.  These  couplings  can  be  calculated  from  (15)  (in  the 
absence  of  a  cover  plane);  the  resulting  nearest-neighbor  terms  in  the  dipole-dipole  coupling 
part  of  the  Coulomb  potential  (17)  are 
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(34) 


In  addition,  there  are  dipole-dipole  couplings  to  neighbors  that  are  further  away  and  that  turn 
out  to  be  larger  than  the  off-diagonal  terms  (/x  ^  v)  in  (34).  The  vibrational  normal-mode  band 
structure  due  to  all  of  these  dipole-dipole  couplings  is  shown  in  figure  4,  with  the  effective 
density  of  states  shown  in  figure  5.  It  consists  of  two  bands,  in  which  neighboring  ions  oscillate 
in-phase  (upper  band)  and  out-of-phase  (lower  band),  and  whose  small  frequency  spread  is 
indicative  of  the  dominance  of  the  nearest-neighbor  coupling  over  all  other  couplings. 

Many  dipole-dipole  couplings  of  the  sort  of  r? rvj  with  /x^v  are  nonzero  in  this 
configuration;  however,  they  do  not  lead  to  effective  spin-spin  couplings  if  the  underlying  trap 
frequencies  along  the  directions  mf  and  m'j  are  strongly  off-resonant  (see  sections  3.3  and  4.2). 
Thus,  neglecting  any  /x  ^  v  couplings,  the  effective  spin-spin  Hamiltonian  that  is  constructed 
from  /x  =  X  Coulomb  interactions  is  approximately 
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(Old)  Aco/co0 

Figure  5.  Density  of  states  of  the  vibrational  bands  (figure  4)  due  to  the 
dipole-dipole  interactions  of  vibrations  along  the  sets  of  axes  in  (33),  shown 
in  figure  3.  Left:  the  six  bands  consisting  of  three  off-resonant  doublets  (//  = 
X,  Y,  Z)  with  central  frequencies  (bare  trap  eigenfrequencies)  split  by  the  golden 
ratio  (see  section  4.2);  do  =  (ci)xthy^z)1/3  and  (jOqy/co  =  0.02  (much  larger  than 
in  a  realistic  experiment).  Right:  zoom  of  one  of  the  doublets.  The  two  bands 
detailed  in  figure  4  are  clearly  visible,  separated  by  ~  4<z>o,  and  show  the  extent 
to  which  the  couplings  of  figure  3  are  dominated  by  the  desired  nearest-neighbor 
couplings.  The  scale  of  the  bands  is  co0 M  =  Q2 /(S7te0do^Md3). 


where  the  first  sum  contains  couplings  between  the  different  lattices  while  the  second  sum 
contains  couplings  within  the  lattices;  numerical  prefactors  for  the  perturbing  terms  are 
used  for  brevity,  as  in  figure  3.  Here,  we  have  assumed  for  simplicity  that  all  ions  are 
simultaneously  pushed  by  the  same  state-dependent  force  with  equal  phase.  The  Hamiltonians 
%Y  and  T-Lz  are  found  from  (36)  through  rotations  by  ±120  ,  and  the  total  effective  spin 
Hamiltonian  is 

EL'  =  —Jxdix  —  JyELy  —  JzELz,  (36) 

where  Jx,  Jy  and  Jz  are  effective  coupling  constants  containing  the  diagonal  coupling  strength, 
the  physical  prefactors,  as  well  as  the  mechanisms  used  for  achieving  these  effective  spin-spin 
couplings  (see  section  4.4).  The  topic  of  whether  or  not  this  slightly  perturbed  Hamiltonian  (36) 
exhibits  the  same  interesting  topological  phases  as  the  ideal  Hamiltonian  (31),  at  zero  or 
finite  [41,  42]  temperature,  is  beyond  the  scope  of  this  paper.  We  mention,  however,  that  if 
the  perturbative  terms  of  (36)  will  be  deemed  too  strong,  they  can  be  reduced  further  by  driving 
the  different  wires  with  different  relative  phases  or  amplitudes  (see  (30)). 

The  presented  configuration  of  trapping  height  and  vibrational  axes  nearly  maximizes 
the  desired  dipole-dipole  couplings  at  the  same  time  as  it  nearly  mimimizes  all  undesired 
couplings.  By  numerical  optimization  we  can  identify  a  configuration  that  performs  a  few 
per  cent  better  than  (33),  but  we  have  not  been  able  to  obtain  an  analytical  description  of  this 
configuration. 

4.2.  SE  trap  design 

To  have  maximally  incommensurate  vibrational  frequencies  along  the  normal-mode  axes  (33), 
we  choose  them  in  the  golden  ratio  cox  :  coY  :  coz  =  4>~{  :  1  :  <p  with  0  =  ( 1  +  V5)/2.  We  use  the 
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Figure  6.  Left:  optimized  rf  (blue)  and  dc  (white)  electrodes  for  the  constraints 
of  section  4.2,  with  no  spurious  traps.  Dimensionless  trap  curvatures  are  k  = 
0.080.  Right:  optimized  electrodes  for  a  honeycomb  lattice  with  out-of-plane 
quadrupole  confinement,  trapping  height  h  =  <7/2  and  cover  plane  at  H  =  50 d. 
Dimensionless  trap  curvatures  are  k  =  0.102.  Coordinates  are  as  in  figure  3. 

algorithm  of  [36]  to  find  an  rf  SE  pattern  that  will  generate  an  infinite  honeycomb  lattice  of 
exactly  such  micro  traps,  with  the  following  constraints: 

•  The  unit  cell  of  the  electrode  pattern  is  defined  by  the  vectors  a  =  d{s/ 3,0,0}  and 
b  =  d{V3/2,  3/2,0}. 

•  The  ion  positions  within  the  unit  cell  define  the  sublattices  LCl  and  L.:  R0o  =  <i{0,  0,  1/2} 
and  R0.=d{V. 3,  1,  1/2}. 

•  The  gradient  of  the  rf  electric  potential  generated  by  the  SEs  must  vanish  at  the  ion 
positions  in  order  to  have  minima  of  the  rf  pseudo-potential. 

•  The  principal  axes  of  the  second  derivative  tensors  of  the  rf  electric  potential  at  the  ion 
positions  are  aligned  with  the  directions  given  in  (33),  with  eigenvalues  proportional  to 
[0_1,  1,  — </>}  in  the  mf .,  mY0%  and  m z #  directions,  respectively. 

•  A  cover  plane  is  located  at  a  height  H  =  50 d. 

The  resulting  electrode  pattern  is  shown  in  the  left  panel  of  figure  6.  It  generates  microtraps  at 
the  desired  positions  with  dimensionless  curvatures  [36]  k  =  0.080  and  no  spurious  additional 
microtraps.  This  is  to  be  compared  with  a  simple  out-of-plane  quadrupole  honeycomb  lattice 
geometry  (k  =  0.102)  as  in  [36]  (see  figure  6,  right  panel),  which  can  potentially  be  deformed 
during  the  experiment  via  dc  electrode  potentials  into  satisfying  the  above  constraints.  Such  dc 
electrodes  might  be  necessary  in  any  experimental  implementation  in  order  to  null  micromotion 
of  the  ions  [43]  induced  by  manufacturing  inaccuracies  and  stray  charges. 

4.3.  Trap  depth  and  trap  loading 

The  depth  of  the  microtrap  lattices  generated  by  the  electrodes  shown  in  figure  6  is  rather 
shallow.  For  the  honeycomb  lattice,  which  is  more  easily  analyzed  due  to  its  p6m  symmetry, 
figure  7  (solid  line)  shows  the  ponderomotive  pseudo-potential  along  a  vertical  axis  through 
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Figure  7.  Total  potential  on  a  vertical  axis  through  any  microtrap  formed  by 
the  electrode  of  figure  6  (right  panel).  The  ponderomotive  pseudo-potential  is 
drawn  with  a  solid  line,  in  units  of  Epp  =  Q2U^/(4MQ,^fd2).  Two  levels  of  dc 
biasing  (either  through  the  rf  electrode  or  applying  the  same  bias  potential  to 
the  dc  electrodes  and  the  cover  plane;  vertically  offset  to  leave  the  trapping 
minimum  unchanged)  are  shown,  in  units  of  Vpp  =  Ew/Q.  For  a  bias  voltage 
V  >  0.1 25  Vpp  the  hexagonal  lattice  microtraps  are  the  only  attractors  (trapping 
zones)  for  cooled  ions. 

any  microtrap,  in  units  of  Evp  =  Q2U2l/(4MQ2ld2).  For  9Be+  ions  (Q  —  +e,  M  =  9  u)  trapped 
with  UT f  =  50  V  and  Qrf  =  2n  x  200  MHz  in  a  lattice  of  d  =  30  /xm  ( h  =d/ 2  =  15  /x m  and 
H  =  50 d  =  1.5  mm),  we  have  Epp  =  4.7  eV  =  5.5  x  104kB  K  and  thus  a  trap  depth  (ion-loss 
barrier)  of  0.001  85Epp  =  8.7 meV  =  101kB  K.  This  small  adiabatic  trap  depth  is  likely  further 
reduced  by  the  breakdown  of  the  pseudo-potential  approximation  near  the  trap  barrier.  In  order 
to  reliably  load  these  microtraps,  we  can  make  use  of  the  cover  plane  (see  section  2.1):  applying 
a  positive  dc  bias  potential  to  the  dc  electrodes  and  the  same  potential  to  the  cover  plane  adds 
an  electrostatic  potential  that  pushes  the  ions  towards  the  rf  electrode  and  increases  the  trapping 
well  depth  (see  figure  7).  Since  this  dc  potential  is  equivalent  to  applying  a  negative  dc  bias 
potential  to  the  rf  electrode,  its  dc  electric  field  at  the  ion  trap  sites  vanishes  (by  construction  of 
the  rf  electrode  shape)  and  it  thus  does  not  induce  micromotion  [43].  We  find  that  applying  a 
small  dc  bias  potential  of  at  least  0.125£pP/ Q  ~  0.6  V  is  sufficient  to  make  the  desired  lattice 
of  microtraps  the  only  minima  of  the  total  potential  (dashed  line  in  figure  7).  By  applying  a 
stronger  bias  voltage,  the  resulting  total  potential  (dotted  line  in  figure  7)  is  deep  enough  to  trap 
ions  produced  by  photoionization  directly  from  a  hot  atomic  beam.  This  bias  will  simultaneously 
cause  the  traps  to  be  shallower  in  the  xy-plane. 

4.4.  Wires  for  magnetic  interaction 

As  described  in  section  3,  effective  spin-spin  interactions  between  ions  require  internal-state- 
dependent  forces  to  be  applied  to  the  ions.  In  the  present  model,  we  propose  to  embed  parallel 
wires  below  the  electrode  plane,  which  generate  local  magnetic  field  gradients  at  the  positions 
of  the  ions  (see  (20)).  A  relatively  simple  periodic  grid  of  two  different  types  of  wires,  indicated 
in  red  and  blue  in  figure  3,  suffices  to  implement  the  spin-spin  interactions  along  all  three 
bond  types  of  the  Kitaev  model.  As  explained  in  section  3,  one  can  induce  pairwise  d{f  ] d{'] 
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and  d'}<°)d'y*,  interactions  by  currents  at  frequencies  that  are  near-resonant  to  a)^ ;  ±  ojx  and 
(o^i  ±  ojy  (M0lmer-S0rensen-type  interactions  [2]);  the  d^]  d(*]  interactions  can  be  driven  with 
currents  that  are  near-resonant  to  a>z  (phase-gate-type  interactions  [5]).  Because  the  three  band- 
manifolds  are  well  separated  in  frequency  (figure  5),  the  dynamics  of  the  three  bond  types  can 
be  driven  simultaneously  by  currents  at  separate  frequencies  that  are  mutually  off-resonant. 

The  geometry  of  the  wires  is  determined  by  the  condition  that  we  need  to  suppress 
the  magnetic  field  at  the  position  of  all  ions,  in  order  to  have  negligible  carrier  interactions 
cp  in  (20),  while  maintaining  useful  field  gradients  that  couple  to  their  target  vibrational 
directions  (33)  to  drive  spin-spin  interactions  on  all  three  bonds  simultaneously  [21,  22].  The 
magnetic  field  of  two  infinite  sets  of  infinitely  long  wires  parallel  to  the  y-axis  as  in  figure  3, 
with  distance  dw  =  dy/3/2  between  wires  of  equal  color,  is 
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where  /x0  is  the  magnetic  constant.  The  field  over  the  blue  wires  (i.e.  at  the  ion  positions) 
vanishes  at  height  hw  if  the  ratio  of  currents  is 
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With  this  current  ratio  the  magnetic-field  gradient  at  the  ion  positions  is 
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(39) 


for  any  n  e  Z.  As  this  magnetic-field  gradient  decreases  rapidly  with  increasing  distance  h w  to 
the  ions,  one  should  place  the  wires  as  close  as  possible  to  the  ions.  On  the  other  hand,  they 
should  not  interfere  with  the  trap  electrodes.  As  a  reasonable  compromise  for  the  following 
estimates,  we  can  assume  that  the  wires  are  below  the  electrodes  such  that  h  w  =  dw  =  hy/3.  The 
actual  hw  in  an  experiment  will  probably  be  dictated  by  constraints  in  the  microfabrication. 

We  choose  the  quantization  axis  of  the  pseudo-spins  of  the  ions  to  coincide  with  its 
associated  bond  direction,  Z  =  A z/d\  however,  any  other  choice  of  Z  will  be  equally  valid,  and 
the  experimenter’s  choice  may  depend  on  the  available  quantization  fields.  For  our  choice,  sz  — 

1 YZ  x  g/rB47r/x3°J’lue  sinh"2  =^-  for  all  ions  on  both  sublattices,  and  with  y^z  =  4n®Mcli 
(the  diagonal  term  of  (34))  the  interaction  strength  (30)  becomes 
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where  g  is  the  effective  g-factor  such  that  the  energy  difference  between  the  |f)  and  f ) 
pseudo-spin  states  in  a  weak  constant  magnetic  field  along  the  quantization  axis  is  A — 


hoon-  gii^Bz.  is  the  current  amplitude  in  the  blue  wires  at  frequency  d)z  +  Sz  with 

\8Z\  <$C  d)z.  With  d  =  30 /xm,  g  =  1,  M  —  9  u,  Q  =  +e  and  coz  =  2jt  x  5  MHz,  this  coupling 
strength  is  Jz  —  7.6  kHz  x  [1^ / A]2\8Z / {2n  kHz)]-2.  To  avoid  sizable  entanglement  of  the 
(pseudo)spins  with  the  ion  motion,  we  need  to  fulfill 
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Jy  |  <  |/j|/2,  in  which  case  the  coupling  constant  of  the 
:  7 2  j2 

■'x  i 


3)  <  I  1/256  (see  [25]  for  discussions  of 


The  hexagonal  Kitaev  model  features  interesting  gapped  phases  with  anyonic  excitations  of 
the  ground  state  for  example  if  |  Jx  \ 
resulting  effective  Hamiltonian  is  7eff  =  7272/(16|/ 
these  phases  and  their  emergence  from  (31)). 

While  the  geometric  prefactor  of  (40)  depends  on  the  details  of  the  model,  its  functional 
dependences  are  expected  to  remain  the  same  for  a  broad  class  of  wire-driven  coupled  pseudo¬ 
spin  models,  in  particular  also  for  Jx  and  JY  of  the  same  system.  The  exact  form  of  the 
interactions  along  the  X  and  Y  bonds  depends  on  the  transition  dipole  matrix  elements  fid  = 
(tl/ill)  which  can  have  components  along  all  spatial  directions.6  The  component  along  the 
quantization  axis  Z  is  relevant  for  tt  transitions  (where  the  Z  component  of  the  total  angular 
momentum  F  of  the  ion  does  not  change,  Am  F  =  0),  while  perpendicular  components  can  be 
used  for  a  transitions  (during  which  the  Z  component  of  F  changes  by  Amp  =  ±1).  The 
coupling  strengths  for  tt  transitions  are  found  by  scaling  (40)  to  be 
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with  currents  of  amplitude  7b(1“j/±ft>x/y')  at  frequencies  07  1  ±  (coX/y  +  8x/y)  that  are  required 
to  drive  Mplmer-Sprensen-type  interactions  [2].  For  a± -transitions  analogous  relations  hold 

A  A  j 

involving  the  projections  of  fid  along  (X  ±  iF) /v2.  We  conclude  that  all  interactions  are  similar 
in  magnitude  and  that  the  current  amplitudes  can  be  used  to  tune  the  effective  coupling  strengths 
Jx,  JY  and  J7  of  the  Kitaev  model  simulator  (36).  To  drive  all  bonds  simultaneously,  a  total 
of  five  alternating  currents  at  different  frequencies  are  necessary;  using  the  largest  allowed 
value  of  1  in  (41),  the  maximum  rms  current  each  of  the  blue  wires  (and  red  wires,  see  (38)) 

has  to  sustain  is  J (7biue(0)  ~  ^5/2  x  0.36  A  x  [5/(27 r  kHz)]3/2.  We  recall,  however,  that  (40) 
and  (42)  depend  very  strongly  on  the  vertical  distance  /1  w ,  and  even  a  small  reduction  in  h  w  can 
substantially  reduce  the  required  currents. 

Expression  (40)  seems  to  suggest  that  for  quantum  simulators  built  with  the  principles 
described  here,  decreasing  the  physical  size  of  the  ion-trap  lattice  (as  given  by  the  length  scale 
d )  will  strongly  increase  the  simulation  speed,  which  is  given  by  the  effective  dynamics  of  the 
particular  quantum  simulator  but  ultimately  proportional  to  the  spin-spin  coupling  strengths. 
However,  a  careful  analysis  of  assumptions  and  constraints,  carried  out  below  for  interactions 
along  Z  but  equally  valid  for  all  other  effective  interactions,  disproves  this  observation.  Firstly, 
if  we  assume  that  avoiding  ion  motion  is  a  significant  experimental  constraint,  a  constant  ratio 
Jz/(h&z)  in  (41)  implies  that  the  effective  spin-spin  coupling  strength  scales  as 


J z  oc 


(43) 


for  a  given  ion  species  and  given  electrode  shapes.  Secondly,  assuming  the  currents  to  be  limited 
by  heat  dissipation,  an  upper  bound  on  must  scale  as  J3/2.  And  lastly,  lower  bounds  for 
the  trap  frequency  can  be  found  in  two  ways:  (i)  to  meet  our  assumption  of  stiff  trapping, 
i.e.  coqz  <$C  d)z,  we  require  coz  ^  d~3/1  x  1 2|/V87rc0M;  and  (ii)  to  use  the  expansion  (26),  we 


6  |  A  A 

In  the  case  of  a  real  spin,  /id—y  g/in(X  —  iF)  (see  the  example  in  section  3.2);  but  what  follows  also  applies  to 
more  general  pseudo-spin  cases  where  we  can  have  fid  ■  Z  ^  0,  for  example  if  the  pseudo-spin  states  are  hyperfine 
states  of  an  ion. 
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require  cooz  <5C  |<5Z|,  which,  combined  with  the  scaling  of  (43)  for  Sz  and  with  the  above  current 
scaling,  implies  a  scaling  of  d~5  for  the  lower  bound  of  0)7.  Together,  these  bounds  imply  that 
the  maximal  achievable  coupling  strength  (43)  scales  only  as  d~ or  even  d2,  depending  on 
which  of  the  frequency  bounds  is  more  stringent.  Since  current  experimental  setups  are  far  from 
reaching  these  lower  bounds  on  &>z,  miniaturization  is  expected  for  now  to  increase  the  coupling 
strength  faster  than  these  estimates,  but  the  optimal  dimension,  where  the  ratio  of  simulation 
speed  and  heating  rate  (anomalous  heating,  scaling  as  d~ 4  [44,  45])  is  maximized,  remains  an 
open  question. 

5.  Conclusions 

We  have  discussed  modifications  to  Coulomb  potentials  and  interactions  of  trapped  ions  due 
to  the  presence  of  trap  electrodes  and  cover  planes.  For  plane  geometries,  we  have  treated 
these  modifications  rigorously,  using  the  method  of  image  charges.  We  have  found  considerable 
deviations  of  the  long-range  behavior  from  that  in  free  space  when  the  relevant  distances  are  of 
the  order  of  ion-to-surface  distances  or  larger.  Moreover,  we  have  developed  a  general  approach 
to  treating  the  effective  spin-spin  interactions  of  ions  trapped  in  a  multi-trap  array  in  the 
stiff-trapping  limit,  where  dipole-dipole  interactions  between  nearest  neighbors  produce  only 
small  corrections  to  the  bare  normal  modes  of  a  given  trap  well.  We  have  shown  that  effective 
coupling  strengths,  and  therefore  simulation  time  scales,  are  determined  by  the  nearest-neighbor 
dipole-dipole  couplings.  As  an  illustration  of  the  versatility  and  power  of  this  stiff-trap-array 
approach,  we  have  discussed  a  quantum  simulation  of  the  hexagonal  Kitaev  model.  We  have 
also  addressed  several  practical  challenges,  including  how  the  trap  depth  of  the  array  may  be 
improved,  so  ions  created  from  a  thermal  source  with  large  kinetic  energies  can  be  trapped. 
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Appendix.  Summary  of  the  used  coordinate  systems 

In  order  to  help  distinguishing  the  various  coordinate  systems  used  in  the  text,  we  summarize 
them  here. 

•  The  laboratory  frame  is  spanned  by  the  unit  vectors 

x  =  {1,0,0},  y  =  {0,1,0},  £  =  {0,0,1}.  (A.l) 

Its  orientation  is  shown  in  figures  1-3.  In  section  2,  laboratory-frame  vectors  are  written  as 
r  =  xx  +  yy  +  zz  with  p  =  sjx2  +  y 2  and  r  =  y/x2  +  y2  +  z2. 

•  The  pseudo-spin  quantization  frame  is  given  by  the  orthonormal  unit  vectors  X,  Y  and  Z, 
where  Z  is  the  quantization  axis.  In  section  4.4,  we  set  Z  =  A z/d. 
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•  The  /  th  ion’s  vibration  around  its  equilibrium  position  is  expressed  in  the  local  coordinate 
frame  mf;  see  (12).  For  the  Kitaev  model,  we  use  the  vectors  given  in  (33):  each  vibrational 
direction  (depending  on  which  sublattice  the  ion  is  located)  is  indexed  by,  and  associated 
with,  one  of  the  spin-space  directions  X,  Y,  Z,  but  this  does  not  mean  that  the  vibrational 
directions  are  parallel  (or  in  any  way  related)  to  the  spin-space  axes. 

•  The  vectors  connecting  neighboring  ions  in  the  Kitaev  honeycomb  lattice  A.x,  Ay  and  Az 
are  of  length  d  and  given  in  (32).  They  all  lie  in  the  plane  of  the  lattice  and  do  not  form  a 
3D  coordinate  system. 
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